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Abstract 



This paper investigates the stabihty of the power-law steady state often observed in 
marine ecosystems. Three dynamical systems are considered, describing the abundance 
of organisms as a function of body mass and time: a "jump-growth" equation, a 
first order approximation which is the widely used McKendrick-von Foerster equation, 
and a second order approximation which is the McKendrick-von Foerster equation 
with a diffusion term. All of these yield a power-law steady state. We derive, for 
the first time, the eigenvalue spectrum for the linearised evolution operator, under 
certain constraints on the parameters. This provides new knowledge of the stability 
properties of the power-law steady state. It is shown analytically that the steady state 
of the McKendrick-von Foerster equation without the diffusion term is always unstable. 
Furthermore, numerical plots show that eigenvalue spectra of the McKendrick-von 
Foerster equation with diffusion give a good approximation to those of the jump-growth 
equation. The steady state is more likely to be stable with a low preferred predator 
: prey mass ratio, a large diet breadth and a high feeding efficiency. The effects of 
demographic stochasticity are also investigated and it is concluded that these are likely 
to be small in real systems. 
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1 Introduction 



It is well established that marine ecosystems often show roughly equal abundances of biomass in log- 
arithmically inc reasing weight i ntervals, when organisms are iden tified by body mass rather than by 



species identity (jSheldon et al 



1972 



Boudreau and Dickid . ll992l ). This is equivalent to a power-law 



for the abundance density as a function of body mass with exponent of app roximately —2. Alter- 
natively, plotting log(abundanc e) against log(mass) gives a "size spectrum" (jSheldon and Parsons . 



1967 



Piatt and Denmanl . Il978l ) which is approximately linear with gradient near to — 1. 



Thi s emp irical pattern has motivated a programme of theoretical research. Silvert and Piatt 
(jl978l : 1 19801 ) developed a size-dependent partial differential equation modelling growth and death 
in a size spectrum, and established the existence of a power-law steady state. The power-law 
steady state has also been shown in s ystern s where predators are allowed to eat any prey smaller 



than themselves (jCamacho and Sola . 



2001 



When predators are assumed to be more selective 
(i.e. eating only a certain range of prey), the existence of a power-law steady state has also been 
proven, using an integro-differential equation for the model instead of a partial differential equation; 
the exponent gen erally depends on assimi lation efficiency, external mortali ty and predator-prey 



2006 



interaction rates (Benoit and Rochet 



Blanchard et al 



20091 : 



Law et al 



20041) . In these and other studies (e.g. lAndersen and Beveij . 



20091 ). the McKendrick-von Foerster equation is commonly 



used. However, a d erivation from a stochastic model of predation leads to a more general equation 
(jPatta et al.l . l2010l ). which we will refer to as the "jump-growth" equation in the following analysis. 
The McKendrick-von Foerster equation is the first order approximation (in an infinite series) to 
the jump- growth equation when prey are typically much smaller than predators. The second order 



appr oximation brings a diffusion term into the McKendrick-von Foerster equation (|Datta et al. 
2010l ). the effects of which have not previously been studied. 



Marine biologists need to understand the resilience of the power-law steady state to perturbations 
caused by fishing and natural phenomena, such as springtime plankton blooms. For instance, it 
has been sh o wn that fishing increases th e temporal variability in abundance of marine species 
( Hsieh et al. . 20061 : Anderson et al. . 20081 ). Fundamental to this understanding are the stability 



properties of the power-law steady state, about which very little is known. We do know from 
recent numerical studies on the jump-growth equation and the McKendrick-von Foerster equation 
that there is a bifurcation from a stable power-law steady stat e to a travelling- wave attractor under 



certain parameter conditions (jLaw et al.l . l2009l : 



Datta et al. 



20 id ) . However , the only s t abilit y 



analysis we are aware of assumed growth to be independent of prey density (jArino et al.l . |2004| ). 
thereby excluding a key predator-prey interaction at the heart of the dynamics. The power-law 
steady state plays a pivotal role in marine ecosystems, and it is essential to understand the factors 
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that contribute to its stability and instablity. 

This paper provides the first detailed stability analysis on the jump-growth equation and its low 
order approximations, the McKendrick-von Foerster equation and the McKendrick-von Foerster 
equation with diffusion. It is also the first analysis of the effects of including the second order 
diffusion term in the McKendrick-von Foerster equation, and of the effects of demographic noise on 
the stable power-law steady state. The results show that the first order approximation is unstable, 
whereas the second order approximation can be stable, and gives a much better approximation 
to the jump-growth equation. The steady state is shown to be more likely to be stable when the 
preferred predator : prey mass ratio is reduced and the diet breadth and the feeding efficiency are 
increased. 

For readers interested in the mathematical derivation of the perturbation equations and eigenvalue 
spectra. Section [2] shows the necessary steps taken. However, for those more interested in the results 
of the stability analyses. Section [3] shows the behaviour of the three models, and reading Section 
12.11 should provide sufficient background reading to understand the different models used. 



2 Analysis of the power-law steady state 



2.1 Three models of predation 



The analysis focuses on perturbations around the power-law steady state of three equations: the 
jump-growth equation ([1]), the McKendrick-von Foerster equation ([2]) and the McKendrick-von 
Foerster equation with diffusion ([3]). These equations describe the rate of change in the density 
of organisms of weight w, which we call (j){w), with dimensions M^-'^L^^, where M is the mass 
dimension and L is the length dimension. This density is with respect to both mass and volume, 
so the number of organisms in a volume V with weigh t between w and w + dw is V(j){w)dw. The 
first equation is based on the jump-growth equation of 



Dattaetal 



mm, 



dcj){w) 
dt 



(^—T{w, 'w')(j){w)<j){w') — T{w', w)4>{w')(j){w) 
+T{w — Kw' , w')(l){w — K'w')(j){w')) dw' — H(j){'w) 



(1) 



T{w,w') is proportional to the feeding rate of individuals of weight w on indivi duals of we i ght w \ 
and < -fC < 1 is the conversion efficiency of biomass from prey to predator (iLaw et al.l . |2009|). 
There are three ways in which a feeding event can result in a change in the density of individuals at 
a given weight w, corresponding to the three terms in the integrand. The first term represents the 
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loss of individuals of weight w due to growth to a larger size (predation of w upon w'), the second 
term the loss of individuals of weight w due to death (predation of w' upon w), and the third term 
the gain of individuals of weight w due to growth from from a smaller size (predation oi w — Kw' 
on w'). Here we have also included a linear natural death rate fi (with the dimension of inverse 
time) to allow for other sources of mortality. 

A Taylor expansion of the third term in the jump-growth e quation in powers o f K gives an infinite 



series of approximations to the full jump-growth equation (jPatta et al.1 . l20ld ). Expanding up to 



and including terms linear in K gives our second model, the McKendrick-von Foerster equation, 

d(j){w) 



dt 



d_ 

dw 



T{w',w)4>{w)(l>{w') dw' (2) 
Kw' T{w,w')(j){w)(l){w') dw' — fj,(p{w) 



and including terms quadratic in K gives our third model, 

d(l){w) 



dt 



-- - J T{w',w)(t){w)^{w') dw' (3) 
d f 

— —— / Kw' T(w,w')(j)(w)(j){w') dw' 
ow J 

+1-^ I {Kw'fT{w,w')4>{w)ct){w') dw'-^i<p{w), 

which we will refer to as the McKendrick-von Foerster equation with diffusion. Note that, as in 
equation ([T]), a linear death rate fi has been included in these two approximations. 

We assume a feeding kernel of the form 

T{w,w') = Aw''s(^) (4) 
\w' / 

where A is the predator search volume per unit mass"" per u nit time, a i s the predator search 



exponent, calculated to have a value of approximately 0.8 (see IWard . Il978l ). and s{w/w') is the 
feeding preference function, centred around some preferred predator : prey mass ratio B. To 
make analytical progress in this paper we assume that a = 7 — 1, where 7 is the exponent of the 
power-law steady state {k. 2). This assumption then has the consequence that the steady state is a 
power-law (see below). In addition, the eigenvalue spectrum can then be written as a closed form 
expression and its properties analysed. Although probably not realistic from a biological point of 
view (discussed in Section H]), the assumption places stability analyses of size spectra on a firm 
mathematical foundation and provides a basis from which exploration of a broader class of systems 
can begin. 
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Section 12.21 defines the power-law steady state for the jump-growth equation ([T]) and its two 
approximations ([2]) and ([3]). Section [2T3l develops equations for the dynamics of small perturbations 
to this steady state and Section 12.41 gives explicit equations for the eigenvalue spectra. In Section 
12.51 the effect of demographic noise on the system at steady state is investigated. Finally, Section 
12.61 incorporates a Gaussian feeding preference for predators. 



2.2 The power law steady state 

The steady state for equations ([T]), ([2]) and ^ is given by 

4){w) = (t)QW~'^ (5) 

where (pQ is a constant. Below, it helps to transform the variable w to a dimensionless log weight 
variable x = ln(it;/it;o) (for some arbitrary weight wq). For analysing the steady state of the jump- 
growth equation ([1]) it is convenient to change the integration variable of each of the three terms 
to the predator : prey mass ratio, which leads to the transformed equation 

dv{x^ ^ f ( 

— A s{e'^) — e°''^v{x)v(x — r) — v{x)v{x + r) 



dt 



_l_g"{»'+V'W)^('2; _ ^(r))t;(x — r — il^ir)) I dr — ^v{x), (6) 



where we have used equation @ for the feeding kernel with a = 7 — 1. Here v{x) has the property 
that e~("^"^^^f (x) dx = (f>{w) dw and has dimensions L~^, A = Awq", and r is the log of the predator 
: prey mass ratio with ^/)(r) = ln(l + Ke~^'). In the transformed jump-growth equation ([6]), the 
steady state is simply given by 

v{x) = vo, (7) 

where vq = (powo^"' is a constant. Substituting this into equation Q we get the steady state 
condition, 

s(e'-) (^_e"'- _ 1 + e"('^+'^('^») dr - t] = (8) 

where r] = ^/(Avq) is dimensionless. This equation implicitly determines the value of the search 
volume exponent a (and thus the steady state exponent 7) for a given choice of the parameters 
K and tj and the feeding kernel s{e^). If we impose the conditions that predators can only feed 
upon prey smaller than themselves and 7^ 0, we can prove analytically that there always exists 
a unique value for a that solves the steady state condition. Without these conditions we verify its 
existence and uniqueness numerically. Setting rj determines the abundance of fish at the steady 
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state, as it contains the constant vq. 
For the McKendrick-von Foerster equation with diffusion ([3]), the steady state condition is 

J s(e^) (^-l + ai^e("-i)^ + a(a-l)^e("-2)'-^ dr-rj = 0, 



(9) 



and for the McKendrick-von Foerster equation without diffusion ([2]), terms of order in equation 
([9]) are ignored. 



2.3 Perturbations around the steady state of the jump-growth 
equation 

We now add a small perturbation to the steady state of the jump- growth equation and observe its 
evolution over time. If the perturbation grows over time, then the steady state is not stable, and 
the system will not stay at the equilibrium; if the perturbation decays, then the steady state is 
locally asymptotically stable. We call the perturbation vq€{x, t) and obtain its evolution equation 
by substituting 

v{x,t) = vo{l + e{x,t)) (10) 

into equation We now assume that we can neglect terms of order because e is taken to 
be very small. For a finite-dimensional dynamical system this can be justified rigorously using 

However in an infinite-dimensional 



Kirchgraberl (119901) 



the Hartman-Grobman theorem (see e.g. 
system this can be more subtle (see e.g. lAulbach and Garavl (|l993l )) and we proceed formally in 
analogy with the finite-dimensional case. We then use condition ([8]) to eliminate terms of order e'', 
so that only terms of remain. This leads to the linearised perturbation equation 



de{x) 



Avo I s(e'^)|^-e°"(e(x) + e(2;-r)) 
-(e(x) + e(x + r)) 



(11) 



We can change integration variables appropriately so that the right hand side of equation (jlip is 
in the form of an integral operator acting on e, 



de{x) 
dt 



Avq I e{m)G{x,m) dm 



(12) 
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where 



Here r = x — m, zi = ln{K/(e^' — 1)), Z2 = ln(e'' — K) and 6 represents the Dirac delta function. 
The integral kernel G{x, m) can be thought of as an infinite-dimensional version of a matrix with 
indices x and m and the task of solving equation (fT2]) thus reduces to finding the 'eigenvectors' and 
'eigenvalues' of this 'matrix'. To define the operator rigorously in the infinite-dimensional case we 
must first restrict the perturbations to the space of square-integrable periodic functions with some 
period L. On this space the operator is compact and thus it is meaningful to speak of its spectrum 
of eigenvalues. In the end we can then take the period L to infinity. 

2.4 Eigenvalue spectra 

We observe that the integral kernel m) depends on x — m only, i.e. it is a convolution kernel. 
Its 'eigenvectors' are given by plane waves, efc(x) = e^^^ , for any /c S M. We refer to k as the 
wavenumber of the plane wave ek{x) and denote the corresponding eigenvalue as X{k). 

The eigenvalues are 

A(A:) = / s(e'') f-e"'' - e'^"" + e"^+("-^'=)V'(0') U + g"*'^'^) dr - rj. (14) 



We refer to the values taken by X{k) as the eigenvalue spectrum. 

A general perturbation can then be expanded in terms of these plane waves and its time evolution 

is 

e{x,t) = y C(A:)e^*^^+^"°^('=)* dA;. (15) 

The expansion coefficient function C{k) is an even function because e(x,t) is real. Notice that if 
any A(fc) has a positive real part then perturbations grow exponentially with time (the factors A 
and vq are positive constants and thus do not affect the coefficient of t), which means that the 
steady state is unstable. 

To derive the eigenvalue spectrum for the McKendrick-von Foerster equation with diffusion from 
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equation (|14p . ip{r) is expanded in powers of K. Taking terms up to and including yields 



m 



s(e'') -e'"'' + K{a - ik)e 



7^2 



{a -ik){a-l- iA:)e(°-2)'^^ (^1 + e"^^'^) dr - r,. 



(16) 

As in equation ([9]), neglecting terms gives the corresponding eigenvalue spectrum for the 
McKendrick-von Foerster equation. 

It is the real part of the eigenvalue that we are interested in, as it is the sign of this that determines 
whether the perturbations grow or die out over time. If, for some wavenumber k Re(A(A;)) is positive, 
then any perturbation containing a component with this wavenumber will grow over time and thus 
the steady state will be unstable. If Re(A(A;)) is negative for all k then all perturbations die out 
over time, and the steady state is stable. 

2.5 Stochastic fluctuations 



The analysis above is concerned with the deterministic jump-growth equation ([T|) and its low-order 
approximations ([2]) and ([3]). In fact, equation ([1]) is t he mean-field equat ion for a stochastic model 
of pairwise encounters between predator and prey (jPatta et al.l . l20ld ). The magnitude of the 
fluctuations due to the demographic noise in the stochastic model is usually a factor of SI 2 smaller 
than the mean-field solution, where is the number of individuals in the system (jvan Kampenl . 
I992I ). For marine ecosystems Q tends to be very larg e, so the fluctuations will b e relatively 
small, but they can nonetheless have important effects (jMcKane and Newmanl . |2005| ). and may 
significantly impact the patterns observed in empirical data. 

In this section we describe how the magnitude of the stochastic fluctuations, and the correlations 
between the fluctuations at different body sizes, can be predicted. In order to make the following 
statements rigorous, one would work in terms of discrete body size intervals, but we work in 
the continuum formally for convenience, which gives the same results. We let n{x, t) be a random 
variable corresponding to the density of individuals of size w = wpe^ at time t . The random variable 
is described by t he stochastic proces s given in previous work (jPatta et al.l . 120101 ). Following the 
method used bv Ivan Kampen (1992), we separate n{x,t) into a deterministic component v{x,t), 
which satisfies the mean-field equations studied above, and a random fluctuation component ^(x, t): 



n{x, t) = ye~"^ ( v{x, t) + n~2VQ^{x, t) 



(17) 



Since the focus of this paper is the stability of the steady state, we restrict attention to the case 
where the deterministic component is at steady state; the results are therefore only relevant in 
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cases where the steady state is stable. The stochastic fluctuations S,{x,t) can be described by a 
Langevin-type equation 



d 



t) = Avo{ / G(x, y)C{y, t)dy + p{x, t) 



(18) 



where the kernel G is given by equation (fT3|) and p{x,t) is a null-mean noise proces s. Details of 
the de rivation of this equation in the general, non-equilibrium setting may be found in 



Datta et al 



(j20ld ). The covariance of noise at two differe nt body sizes is d escribed by a covariance kernel 



B{x,y) = {p{x,t)p{y,t)) , which is given by (see 



Datta et al 



201 



J: 



B{x,y) = e"^^+y^ I (^f{x,y, 



z) - f{x,z,y) - fiz,y,x) 



+ S{x-y) I ( fix, z, z') + Zifli^^ dz' I dz. 



(19) 



where 



and 



fix, y, z) = e-°(^+^) ikix, y, z) + Hy, x, z)) 



kix,y,z) 



e s 



°^«(e^-^) Siz-x-tl^ix-y)). 



(20) 
(21) 



The covariance (^(x, i)^(y, t)) of the fluctuations at logarithmic body sizes x and y satisfies 

^^{Cix,my,t)) = Avo iGix,zmz,my,t))+Giy,zmz,mx,t)))dz + Bix,y)^ . (22) 

In the steady state, the time derivative on the left-hand side vanishes and thus the covariance 
function {S,ix,t)(,iy,t)) can be calculated by setting the right-hand side to zero which results in a 
linear equation to be solved. We present the numerical results of this in section 13. 7i In order to 
verify these, we also carry out stochastic simulations of the number ni(t) of individuals in log weight 
bracket [xi,Xi+i] for —4 < x < 4 (outside this size range, the spectrum is assumed to remain at 
steady state). We approximate the number Rij of individuals in br acket i that eat an individual in 
bracket j during a short time dt as a Poisson random variable (see 
The mean of Rij is given by 

y"^r(woe^* , woe^^)ni it)nj it)5t. 



Datta et al 



20ld . for details). 



(23) 



The fluctuation ^(xi,t) is computed from the difference between ni(t) and its equilibrium value. 
The covariance (^(xj, t)^(xj, t)) is then obtained by averaging ^(xj,t)^(xj,t) over a large number 
of successive time points. In the long term, this gives the same result as the ensemble average of 
^(xi,t)^(xj,t) provided the stochastic process is ergodic. 
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2.6 Gaussian feeding preference 



Organisms do not eat indiscriminately; here we assume that they feed at some preferred prey size 
(in relation to their own size), and a range of sizes around this preferred size. To reflect this, a 
suita ble preference function is a Gaussian feeding p reference, with peak at j3 and width proportional 
to a (lAndersen and Beyerl . l2006l : iLaw et al, This can be represented by the following form 



for s{e^), 



1 



(24) 



In theory this function allows predators to eat prey larger than themselves (i.e. is non-zero for 
r < 0), although for realistic sets of parameter values s(e'') is typically negligible for r < 0. 

The eigenvalue spectrum for the jump- growth equation (|14|) with the Gaussian preference func- 
tion (|24|) unfortunately does not have a closed form. In contrast, the eigenvalue spectra for the 
McKendrick-von Foerster equation without and with diffusion can be determined analytically. 
Defining 



Rn = (a - n) ( /3 + -(T^(a - n) 



(/3 + o^ia - n)) 



(25) 
(26) 



and taking the steady state condition ([9]) into account, the eigenvalue spectrum for the equation 
with diffusion is 



Re(A(A;)) = e'^' 



cos{kf3) + Ke^^ (acos(/i) - ksm^h)) 



(27) 



2 



+^e^^ ((a(a - 1) - k'^) cos(/2) - k{2a - 1) sin(/2)) 



The diffusion term is removed by excluding terms of order in equation (j27p . An important 
difference between the two approximations is that there must always exist values of k for which 
Re(A(A;)) is positive in the eigenvalue function for the McKendrick-von Foerster equation. Conse- 
quentially the McKendrick-von Foerster equation will never give a stable spectrum. In contrast, 
the McKendrick-von Foerster equation with diffusion contains a non-oscillatory term in k which is 
negative and increases in magnitude as k increases. This has the effect of making the real parts 
of the eigenvalues more negative for higher values of k. For both approximations the oscillatory 
terms are damped by a factor of e~2'^ . Equation ([27|l is analysed in greater detail in Section [3] 
to explain observed patterns in the behaviour of eigenvalue spectra when altering parameters. 



10 



0.8 p 
0.6 
0.4- 
0.2 
^ - 
1 -0.2 - 
-0.4- 
-0.6 
-0.8 



(a) 




Figure 1: The eigenvalue spectra for the (a) jump-growth equation (JGE), (b) McKendrick- 
von Foerster equation (MvF) and (c) McKendrick-von Foerster equation with diffusion (MvF- 
D) when using a Gaussian feeding preference. Note that t] has been set to give a steady state 
exponent of roughly 2.3 for all three spectra. Parameter values K = 0.2, /3 = 5, cr = 1.5, r] = 
0.290, 7 = 2.30. 



Using the steady state condition ([8|), it can be shown for the jump-growth equation that Re(A(0)) = 
r/. This result also applies to both of the approximations. Thus, for any positive r/, Re(A(fc)) must 
be positive at A; = 0, and as A(fc) given in equation (jl4p is continuous, there exists a neighbour- 
hood around k = where Re(A(A;)) > 0. Therefore there will be a range of wavenumbers k for 
which perturbations e^^^ will destabilise the steady state. However, we only expect pur model to 



be re alistic for a range of body weights spanning around 12 orders of magnitude (I Cohen et al 



20031 ) and therefore should ignore perturbations with a wavelength longer than this, i.e. those with 
wavenumbers smaller than about k ~ 0.2. 



3 Results 



3.1 Eigenvalue spectra of the three models 



To evaluate the eigenvalue spectra, we use the Gaussian feeding preference (iMjl , for given values for 
the parameters K, /3, cr, ry. Where p ossible we keep these param eters biologically reasonable and 
close to values from previous studies (jAndersen and Beveri . l2006l ) . We use values oi K = 0.2, /3 = 5 
and (T = 1.5 as a base parameter set, and investigate the effects of changing these parameters. 
For this base parameter set, the steady state exponent 7 is equal to 2.27 when r] = (i.e. no 
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external mortality) and the value of 7 increases with r]. The values of 7 used i n the numerical plots 



mostly lie in the empirical range of 2.2 to 3.25 reported bv lBlan chard et al.l (|2009l ). Values of the 
wavenumber k are taken over a range from to 30, as the interesting behaviour of the eigenvalue 
spectra is seen in this frequency range. Note that the expressions for Re(A(/c)) are even in k for the 
three models, so the plotting of negative values of k is unnecessary. We often plot the eigenvalue 
spectra over a logarithmic A:-axis to make it easier to see the details at small k. 

Examples of the eigenvalue spectra of the jump-growth equation and its two approximations (all 
computed numerically using the preference function (j24p ) are compared in Figure [TJ All three 
spectra are close to rj for small k, as expected from Section 12.61 Both approximations are close 
to the jump-growth equation for low values of k, but as k gets larger only the McKendrick-von 
Foerster equation with diffusion follows the jump-growth equation closely. This is expected from 
equation ([27|) because the diffusion term is needed to make the eigenvalue spectrum more negative 
with increasing k. Adding the diffusion term gives a better approximation to the full jump-growth 
model. Thus the properties of equation (\27\i will be used to gain insight into the behaviour of the 
eigenvalue spectra of the jump-growth equation in the subsequent sections. 

The power-law steady state is unstable for all three models in this example, because all three 
spectra contain eigenvalues with a positive real part (the maximum occurring at /c ~ 0.861). The 
tendency for unstable steady states to emerge in our analysis will be discussed in Section [H 

The different behaviours of the two approximations are not just limited to a Gaussian feeding 
preference; similar results have also been obtained when using a step function for the feeding 
preference. This has the form 

= if/'--<'-</' + - (28) 

I otherwise 

and is a rectangular kernel, with midpoint /3, width 2a and height 1/(20"). It is worth noting that 
although behaviour similar to Figure [T] is observed, the oscillations are not damped exponentially, 
and oscillations are observed at all values of k. 



3.2 Stable and unstable steady states 

For some sets of parameter values, the steady state is stable. Figure [2] gives an example, obtained 
by allowing a low preferred predator : prey mass ratio /3, a high efficiency K and a relatively large 
diet breadth a. This example is chosen to illustrate the point that the eigenvalue spectrum for 
the McKendrick-von Foerster equation can be misleading; the spectrum for the McKendrick-von 
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Figure 2: Eigenvalue spectra for the McKendrick-von Foerster equation and McKendrick-von 
Foerster equation with diffusion, compared to that of the jump-growth equation. Parameter 
values K = 0.8, /3 = 1, a = 0.35, 77 = 0, 7 = 2.11. 

Foerster equation without diffusion peaks at 1.19, whereas for the equation with diffusion and the 
jump-growth equation Re(A(A;)) < for all k. The spectrum is stabilised by the non-oscillatory 
term introduced by the inclusion of terms of order K^. The diffusion term contributes to stability 
and the effect of this is great enough to make a qualitative difference to the calculated stability 
of the steady state. As predicted in Section 12.61 the McKendrick-von Foerster equation gives an 
unstable spectrum for any choice of parameter values. 



3.3 Time evolution of perturbations 

To show the consequences of stable and unstable steady states on the dynamics, we can examine 
the behaviour of a local perturbation to the size spectrum, and observe its time evolution. Assume 
a Gaussian perturbation with initial form 

e(x,0) = z^e ^, (29) 
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Figure 3: The time evolution of a Gaussian perturbation, q = 0.2, when using (a) a stable 
eigenvalue spectrum (parameter values K = 0.8, /3 = 1, a = 0.35, r/ = 0, 7 = 2.11), and 
(b) an unstable eigenvalue spectrum (parameter values K = 0.2, /? = 5, a = 1.5, r] = 0.290, 
7 = 2.30). 



where is a smah constant and ? dictates what range of body sizes in the size spectrum are effected 
by the initial perturbation. This can be expanded in plane waves, rewriting €{x, 0) as 

/oo 
6-2^''=%*'=^ dfc. (30) 

where u = (i^?) / ^/2tt. The time dependence of this perturbation then has the fohowing form: 

/oo 
-00 

We set ? so that the perturbation covers about one size unit on the x-scale, and choose units so that 
Avq = 1. We choose to centre our perturbation around x = without loss of generality. Plotting 
the time evolution both for a stable spectrum (Figure E]) and an unstable spectrum (Figure [T]) using 
the jump-growth equation gives the two behaviours shown in Figure [3l 

For both plots, the initial perturbation moves along the x-axis over time, as the organisms it 
contains feed on smaller organisms and grow. In the case of a stable spectrum, the perturbation 
gives rise to smaller peaks either side of the initial perturbation, and these all die out over time. 
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Figure 4: The eigenvalue equations for the jump-growth equation with varying logarithm of 
the preferred predator : prey mass ratio (3. Parameter values a = 1.5, K = 0.2, i] = 0. 

tending to zero across the whole range of x. In the case of an unstable spectrum, the peaks grow 
over time. They develop into waves with wavenumber k, where k is the most unstable node of 
the eigenvalue spectrum. Thus, in the case of Figure [3)3, where k = 0.861, the wavelength of the 
peaks is seen to be around {2Ti)/k. Over time the peaks grow in magnitude but maintain their 
wavelength. The speed at which the perturbation moves through the size spectrum is determined 
by Im(A(A:)). 



3.4 Changing the preferred predator : prey mass ratio 

Figure U] shows the effect of increasing the logarithm of the preferred predator : prey mass ratio /3 
on the stability of the jump-growth equation. The maximum real part of the eigenvalues increases 
as /3 increases, the steady state going from stability when (3 = 1 (i.e. Re(A(A:)) < for all k), 
to instability for the larger values of /3. This is in keeping with previous numerical results, where 
i ncreasing led to a bifurcation from the power-law steady state to a travelling wave attractor 



(|Law et al.l . l2009l ). although the two results should not be directly compared because in earlier 
work the assumption a = 7 — 1 was not imposed. The changes seen in Figure H] as /3 is increased 
can be understood in terms of equation ([27l) . where f3 occurs both in the Rn exponential terms and 
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Figure 5: The eigenvalue equations for the jump-growth equation, with varying feeding 
efficiency K\ k denotes the location of the most unstable node of the spectrum. Parameter 
values 13 = 5, a = 1.5, 77 = 0.290. 

in the I„ cosine and sine terms. In (3 acts to dampen the waves more as it increases, and in /„, 
j3 acts to reduce the period of the waves as it increases. Both these changes are visible in the figure. 
Some decrease in the exponent of the power-law steady state is also evident with increasing (3 in 
Figure m We interpret this in biological terms as an outcome of less biomass being lost from the size 
spectrum as (3 increases, because biomass is inefficiently consumed fewer times during its passage 
along the spectrum. Note that a has been held constant this figure, so that as j3 is increased, the 
mean of the predator : prey feeding distribution increases but the variance remains constant. 



3.5 Changing the feeding efficiency 

Figure [5] shows the effect of changing the feeding efficiency K on the eigenvalue spectrum. To 
understand Figure [H it helps to consider the limiting case of -fC — )• 0. Although unrealistic, because 
it implies no growth of organisms, the eigenvalue spectrum in equation (I27p is then simply a damped 
cosine wave: Re(A(A;)) = e~2'^ ' cos(A;/3). Consequently, the most unstable node k must be the 
first peak of this wave, which occurs at A: = 7r//3, equivalent to k = 0.628 with the parameter values 
in Figure [5l We observe in Figure [5] that, for small K {1 x 10^^), the value of k (0.655) is close to 
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Figure 6: The eigenvalue spectra for the jump-growth equation with varying diet breadth: (a) 
a = 0.25 {k* = 9.60), (b) a = 0.5 {k* = 6.02), and (c) a = 1 {k* = 3.41), where k* denotes the 
largest value of k for which Re(A(A:)) > 0. Parameter values K = 0.2, /3 = 5, 77 = 0, 7 = 2.27. 



this limiting value. 

Corresponding to the node at k = tt/P, there is a dominant eigenfunction with a wavelength 2/3. 
This can be understood in biological terms as a straightforward consequence of the predator-prey 
interaction. A pulse perturbation from steady state that increases the density of predators at some 
size lowers the density of prey times smaller than themselves. This in turn reduces the mortality 
rate on the prey's prey e^^ times smaller than the predators, allowing their density to increase. 
This leads to the wavelength 2/3. 

Figure [5] also shows that, as K increases, k grows and Iie{X{k)) gets smaller. In other words, 
perturbations from the steady state grow more slowly and have wavelengths less than 2/3 as K 
increases. In this case, a pulse increase in predator density at some body size does not remain at 
the same position in the size spectrum as time goes on. The predators grow as they eat, and their 
preferred prey body size moves along with them. This mitigates to some extent the destabilizing 
feedback of slow (or absent) predator growth that would continue to reduce the density of prey 
approxima t ely times smaller than the predator. These results help explain the observation of 



Lawetal 



((20091) that perturbations tend to have a wavelength less than 2/3. Notice also that the 
exponent of the power-law steady state becomes substantially smaller as K increases, because more 
biomass passes along the size spectrum to large organisms. 
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3.6 Changing diet breadth 



It has been shown in earUer numerical studies that, by making the diet breadth more narrow 
(i.e. decreasing a), the power-law steady state can be come unstable, leading to travelling waves of 



Datta et al 



2010|). In the 



abundance that move along the spectrum with time (|Law et al.l . l2009l : 
extreme case of a feeding kernel where predators only eat prey of the exact preferred mass ratio, 
and of no other weight (using a Dirac delta function of the form s{e^) = 6{r — /3) as the feeding 
preference), the steady state can be shown always to be unstable (proof not given here). 

In Figure[6]we investigate the effect of increasing the diet breath a. As a increases, the amplitude 
of oscillations at low values of k decreases, and the range for which Re(A(A;)) has positive values 
becomes narrower; the largest value of k for which Re(A(A;)) > /c* is seen to decrease as a increases. 
This is consistent with equation ([27|) . because increasing a will cause the oscillations to be damped 



sooner by the e 



term. Note that it is the change in a which is causing the change in the 



spectrum and not the steady state exponent; 7 remains at a value of approximately 2.27 in each 
case. 



3.7 The effects of demographic stochasticity 

As explained in section [23} we can calculate the equal-time covariance function {^{x,t)(^{y,t)) for 
the fluctuations in the steady state. This function describes how the fluctuations due to demo- 
graphic stochasticity are correlated at different weights at steady state. It is obtained by solving 
the linear integral equation obtained by setting the time derivative to zero in equation (j22p . To 
perform the calculation we used discrete weight brackets, so that the integral equation becomes a 
matrix equation, which we solved numerically. 

Using the same parameter values as in Figure [2l the covariance function calculated from (|22p 
(solid curve) and from stochastic simulations (points) is plotted in Figure [71 The simulations used 
a system size of = 10^ and were restricted to the size range —4 < x < 4; simulating a wider size 
range for sufficient time was computationally prohibitive. The simulation data become increasingly 
noisy as x increases, due to the decreasing number of individuals in a weight bracket. Nevertheless, 
the simulation results show good agreement with the solution of equation (|22|) . 

The graph in Figure [7] decays exponentially with distance, a typical feature of covariance func- 
tions. Superimposed on the decay is an oscillation with a wavelength of approximately 2/3, generated 
by the non-local predator-prey interaction. The reason for the oscillation is that a positive fluctua- 
tion at X = gives more food, faster growth and a negative fluctuation near /3, which in turn gives 
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Figure 7: The covariance function {^{0,t)^{x,t)) for the fluctuations around a stable steady 
state due to demographic stochasticity (sohd curve - solution of equation ( 12^ : points - 
stochastic simulations). This measures the covariance between stochastic fluctuations at 
log body sizes and x; the covariance between fluctuations at log body sizes x and y, 
{^{x,t)^{y,t)), is equal to e"^(^(0, — x,t)). Parameter values K = 0.8, (3 = 1, a = 
0.35, = 0, 7 = 2.11. 



less food, slower growth and a positive fluctuation near 2/3. 



4 Discussion 



We have presented a local stability analysis of the power-law steady state of marine size spectra. 
The approach has s ome rese i nblance to the local stabilit y analyses of steady-state food webs widely 
applied in ecology ([Murravl . |2002| : iRoonev et al.l . |2006| ). However, instead of having nodes repre- 
senting a finite number of species, the analysis here uses a continuous weight range corresponding 



to an infinite number of "nodes", and this gives a continuous sp ectrum of eigenva. 



ization of the eigenvalue spectrum has been carried out before (jArino et al.l . |200J); the difference 



ues. Character- 



here is that we explicitly link growth of the organisms to predation, which we think is a useful step 
towards reality. 

To do the analysis, the predator search exponent a and steady state exponent 7 have been 
set so that q = 7 — 1. In addition, we assume that the rate for predation-independent death is 
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independent of body weight. These assumptions imply that the dynamics of small perturbations are 
described by the convolution operator given in equation (fT2]) . leading to a simple time dependence 
of the perturbations in terms of an expansion in plane waves, given in equation (jlSp . In general 
this assumption would not be appropriate in ecological communities. The reason for using it here 
is that we believe it is valuable to have analytical results for this special case before beginning 
numerical explorations of conditions closer to those in nature. 

The benchmark for the analysis is a jump-growth e quation, obtained as the large-system limit 
of an underlying stochastic predation-growth process (jPatta et al.l . 120101 ). Importantly, the eigen- 
value spectrum of the well-known, first-order approximation, the McK endrick-von Foerster equation 



( Andersen and Bever 



2006 



Law et al. 



200S 



Blanchard et al 



20091 ). exhibits a systematic depar- 



ture from that of the jump-growth equation: the real parts of the eigenvalues of the former tend 
to zero as wavenumber increases, whereas those of the latter become increasingly negative. There- 
fore in our analysis the eigenvalue spectrum of the McKendrick-von Foerster equation must always 
contain eigenvalues with positive real parts, and must always have an unstable steady state. 

In contrast to the first-order approximation, the eigenvalue spectrum of the second-order approx- 
imation, obtained by adding a diffusion term to the McKendrick-von Foerster equation, contains a 
negative term that is quadratic in the wavenumber, which makes the real parts of the eigenvalues 
much closer to those of the jump-growth equation. The diffusion term is potentially important. One 
consequence of it is that there can be eigenvalue spectra for which Re(A(/c)) < for all wavenum- 
bers A; > 0, implying local stability of the steady state. This is with the caveat that the eigenvalue 
spectrum tends to the natural death rate t] as the wavenumber tends to zero, so perturbations with 
sufficiently low wavenumbers (long wavelengths) could still destabilize the steady state. 

The second-order approximation with diffusion has not previously been used, but would be 
worth considering i n the f uture when the full jump-growth equation cannot be used. Interestingly, 
Benoit and Rochetl (j2004l ) found they had to include a diffusion term in numerical integrations 
of the McKendrick-von Foerster equation to obtain a solution in the absence of natural mortality, 
although they stated that they did not understand why this should be so. How serious the omission 
of the diffusion term is in practice depends on the wavenumber k at which the eigenvalue spectrum 
peaks, because it is this wavenumber that dominates the solution in the long term. If the peak 
occurs at sufficiently small k, the effect of the negative second-order term in equation (j27p is small, 
and the standard McKendrick-von Foerster equation is reliable (Figure [1]) . If the peak occurs at 
large k, the negative second-order term in equation (j27p becomes significant, and inferences about 
stability from McKendrick-von Foerster equation may not be reliable (Figure [2]). The second-order 
equation with diffusion itself becomes a poor approximation if the feeding preference function is 
set such that predators are often smaller than their prey, because the Taylor expansion of the 
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jump-growth equation on which it is based is no longer convergent (jPatta et al.1 . l20ld ). However, 
in reahty predators are almost always larger than prey, so this is not likely to be an issue. 



Key parameters for locating the peak of the eigenvalue spectrum with respect to k are the 
logarithm of the preferred predator : prey mass ratio f3, the efficiency of mass transfer from prey 
to predator K and the diet breadth a. The results in Section 13.51 suggest that predator-prey 
interactions would typically restrict the wavenumber k at the peak to be greater than vr//3. Overall, 
to get the peak of the eigenvalue spectrum at a low wavenumber where the McKendrick-von Foerster 
equation works best, Ke~^ must be small, i.e. growth increments of predators must be small. As /3 
is made smaller and K is made larger, the McKendrick-von Foerster approximation works less well, 
because it misses the stabilizing effect of the diffusion term. The diet breadth a, also affects the 
shape of the eigenvalue spectrum, the main effect in equation (j27p being to dampen the oscillations 
in the real parts of the eigenvalues (Figure [6]). In so doing a has the potential to shift positive 
peaks below Re(A(A;)) = 0, and hence to change an unstable steady state into a stable one. This is 
consistent with the results of earlier studies which have shown the stabilizing effects of broad diets 
20091 : 



( Law et al. 



Datta et al. 



20101). 



Random variability from one individual to another in, for example, the number and size of 
prey items encountered over a period of time, can have important effects in systems such as the 
one studied in this paper. This intrinsic demographic stochasticity is distinct from environmental 
stochasticity, which is not included in the current model. Usually, the relative magnitude of fluc- 
tuations due to demograp hic stochastic i ty is p roportional to (7^^/^, where is the total number of 
individuals in the system (jvan Kampenl . ll992l ). However, interaction between the natural frequency 
of the mean- field system and intrinsic varia bility, which acts at all frequen cies, can cause resonant 
amplification of demographic stochasticity (jMcKane aiid Newmanl. l2005l). S tochastic effects can 



also cause switching between diffe rent solutions ( see e. g. ISamoilov et al.l . l2005l ). a phenomenon that 
cannot be investigated using the Ivan KampenI (jl992l ) approach taken in this paper, which deals 
with fiuctuations about a mean- field solution. As seen in Figure [71 correlations between simultane- 
ous demographic fluctuations at different body sizes do exist. These arise from the predator-prey 
interactions described by the feeding kernel. An investigation of correlation of fluctuations across 
different times is beyond the scope of this paper. However, neither resonant amplification nor 
switching were observed in stochastic simulations, even for relatively small system size ri, and the 
simultaneous correlations observed are likely to be very small for realistic system sizes. 

A feature of the stability analysis here is that the parameter values required to achieve stability 
are outside the range likely to apply in marine systems (e.g. K = 0.8 in Figure[2]). As stated above, 
earlier numerical integrations using the McKendrick-von Foerster equation have led to stable steady 
states using realistic sets of parameter values. There are, however, some important differences 
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between the present analysis and previous work. First, re al size spectra span a finite range of 



body sizes, about twelve orders of magnitude being realistic (jCohen et al.l . l2003l ) . This means that 
perturbations with very long wavelengths cannot occur, and corresponding to this, the wavenumber 
k cannot be less than about 0.2. Second, the finite range calls for lower and upper bounds which 
are not used here. Imposing such bounds removes the exact power-law steady state, and the 
boundary conditions themselves influence the stability of the steady state. Third, the constraint on 
parameter values needed to achieve = 7 — ! may exclude those values likely to lead to stability. 
The present study is best thought of as throwing light on the role that mortality, predation and 
growth play in determining stability of the power-law steady state. Other processes also leave their 
own footprint, and some o f these increase the parameter space in which stable steady states arise 



(jCapitan and Deliud . l20ld ). 



Nonetheless, at a qualitative level, the results here are consistent with earlier observations that the 
steady state of marine size spectra undergoes a bifurcation from stability to instability as predator : 
prey mass ratio is increased and as diet breadth is decreased. The results here indicate that this is a 
Hopf bifurcation as a complex conjugate pair of eigenvalues cross the imaginary axis. Even without 
taking other major life processes into account, the analysis makes clearer what kinds of ecosystems 
are more vulnerable to external disturbances such as those caused by fishing and climate change. 
Further research should expand upon this, to better understand marine ecosystem dynamics, and 
better predict the potential consequences of perturbing seemingly robust ecosystems. 
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